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We present analytical expressions for the Fourier analog of the CMB three-point and four-point 
correlation functions, the spatial bispectrum and trispectrum, of the Ostriker and Vishniac effect in 
the linear and mildly nonlinear regime. Through this systematic study, we illustrate a technique to 
tackle the calculation of such statistics making use of the effects of its small-angle and vector-like 
properties through the Limber approximation. Finally we discuss its configuration dependence and 
detectability in the context of Gaussian theories for the currently favored flat ACDM cosmology. 
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O , I. INTRODUCTION 

(N 

In recent years, with the prospect of the increase in the sensitivity and angular resolution of the forthcoming 
Cosmic Microwave Background (CMB) satellite and interferometry experiments [ESSHj efforts have been driven 
£SJ ■ to the theoretical study of the secondary anisotropies contributions to the temperature fluctuations on arcminutc 
t-H \ scales and below. While the primordial anisotropies from recombination are thought to be well understood secondary 
anisotropies from reionization are not. 

As is well known, the current favored inflationary model of structure formation predicts a nearly Gaussian probability 
distribution for the primordial density fluctuations. In this case, the CMB is completely described by its two-point 
correlation function or power spectrum in Fourier space. All higher-order correlations can be expressed in terms 
of it. Primordial nonlinearities and secondary effects introduce deviations from Gaussianity, producing a detectable 
signal in both the power spectrum and higher-order statistics. Recent work provides the theoretical background for 
t-H | the calculation of estimators of these higher-order statistics 0, |(| and constrains possible non-Gaussian primordial 
contributions to the bispectrum and the trispectrum on degree and sub-degree angular scales using actual data 
[zl IE El OH- The interest is now in forecasting the expected signals on smaller scales due to secondary anisotropies, 
checking whether they are detectable and understanding how they can be separated from each other and from the 
Q^' primary anisotropics in light of the future data. 

The Ostriker and Vishniac (OV) effect was found to be the dominant linear secondary contribution to the CMB 
anisotropies below the Silk-damping scales at the arcminute level ^ is caused by Thomson scattering off of the 
CMB photons by moving electrons during the initial phase of reionization. It has the advantage of taking place during 
d • the linear regime of structure evolution and of being a small-angle effect enabling one to obtain analytical expressions 
for its higher-order correlation functions in the small- angle limit. Because of the highly predictive power of linear 
theory, any measurement of such statistics would be a sensitive probe of the reionization history of the universe, 
difficult to disentangle in measurements from nonlinear contributions. Several derivations for its power spectrum have 
been carried out j 1 2L Il3l 1 1 11 ] but no one has fully addressed the calculation of its bispectrum or trispectrum. 

It is then timely to obtain these expressions and to qualify and quantify them. We therefore extend and detail 
previous techniques used for the calculation of the OV power spectrum to the calculation of its higher-order statistics. 
As it will be shown, for the particular case of fields which are vector-like in nature, such as the OV effect, even moments 
will dominate over odd moments, making the trispectrum a much more sensitive statistics than the bispectrum. 

Given the low redshifts of formation of structure, it is interesting to consider whether nonlinear effects can further 
enhance these statistics. So we will extend our study to the weakly nonlinear regime, allowing us to probe the most 
natural extension of the OV effect to nonlinear scales, the so-called kinetic Sunyaev-Zel'dovich (kSZ) effect. On small 
scales, both arise from the density modulation of the Doppler effect from large-scale bulk flows. 

We review the relevant properties and parameters of the adiabatic cold dark matter (CDM) cosmology for structure 
formation in Sec. Ill Al In Sec. Ill Bl we review the theory of the OV effect and in Sec. Ill CI we discuss the basic 
statistical properties of a general field through its n-point functions. In Sec. Ill Dl we show how the homogeneous 
theory of turbulence combined with the Limber approximation enables one to infer the dominant contribution among 
n-point statistics of a vector field effect like the OV effect. In Sec. Ill El we introduce the standard formalism of the 
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signal-to-noise calculation by means of the Fisher matrix. For general consistency, in Sec. Ill Fl we apply our method 
in detail to the calculation of the OV power spectrum as well as its nonlinear extension. In Sec. 11111 and in Sec. I1VI 
we present the steps of the calculation of the OV bispectrum and of the trispectrum respectively and its nonlinear 
counterparts. We also present the results. Finally in Sec. we conclude. In the Appendix A we generalize the Limber 
approximation to the 3-point and 4-point correlation functions. This may be useful for other cosmological studies. 



II. GENERAL CONSIDERATIONS 



A. Cosmological model 

We work in the context of the adiabatic cold dark matter (CDM) family of models. In units of the critical density, 
fio is the contribution to the non-relativistic-matter density, is the contribution to the baryonic matter density, 
Da is the contribution of the cosmological constant and Ho = 100 ft. km sec -1 Mpc -1 is the Hubble constant today. 
The Friedmann equations for the evolution of the scale factor of the Universe, a(i), are then 



- = h q e(z) = H Q ^n (i + zf + n A + (i - q - n A )(i + z) 2 , (i) 



- = H 2 [n A -n Q (i + zf/2i (2) 

where the over-dot denotes a derivative with respect to time. The scale factor is chosen such that cioHq = 2c. 

Useful measures of distance (and time) are the conformal distance (and conformal time). If an observer is at 
the origin z = then an object at redshift z is at a comoving distance, w{z) — | J* e{z i ) an< ^ a ^ a time t{z) = 
77o (i+z')E(z') • conformal time is then obtained from d-q = dt/a, such that the comoving distance to the 

horizon is the conformal time today, crjo = w(po). 

If the CDM density contrast at comoving position w at time t is S(w, t), then the power spectrum P(k, t) is denned 
by the expectation value over all realizations < d(k,t)6* (k' ,i) >— (2w) 3 Sj D (k — k')P(k,t) where 8^ D is the Dirac delta 
function. In linear theory, 5(w,t) = So(w)D(t)/ D(to), where to is the age of the Universe, 5q(w) = S(w,to), and the 
growth factor, as a function of redshift, is [15| 

The power spectrum is given by P(k,t) = P(k)(D/ Do) 2 , where Do = D(to)- For P(k) = P(k,to) we use 

P(k) = ^4(fc/2)"T 2 (fc p Mpc/ftr), (4) 

where T{q) is the CDM transfer function, k p = kao = kHo/2c is the physical wave number with our conventions and 
T, the shape parameter, is defined as 0] T ~ ^o(/i/0.5) exp(— Qt, — Qb/^o)- As we chose aoH = 2c, there is an extra 
factor of 8 in the denominator in Eq. (0}. For the transfer function, we use the Bardeen et al. 17] fitting formulae for 
CDM models instead of the improved version of Eisenstein and Hu Ha. to facilitate comparison with previous work. 
For Sh, we take the fits to the Cosmic Background Explorer given in [ljj. 

In linear theory, the continuity equation relates the Fourier components of the velocity field and the density field 

m t) = ^ | mi t) = ^ §- km. ( 5 ) 



A useful relation [T5j is 

Dad 5fio a (1 + z) 2 



D a a 2 a[E(z)] 2 D(z 



(6) 



Although we maintain generality in the derivations, we illustrate our results with a flat model, the ACDM model. 
The parameters for this model are f2 = 0.35, = 0.05, Qa — 0.65, h — 0.65 and spectral index n = 1. Concerning 
the reionization contribution we consider two reionization histories, both assuming steep reionization with ionization 



3 



fraction x e = 1 and Az r /(1 + z r ) = 0.1. In the first one, reionization takes place at z r = 8. The second one assumes 
z r = 17 and relics on the latest results from the Wilkinson Microwave Anisotropy Probe (WMAP) experiment (see 
below). Note that in an open or closed universe one replaces in the factors of r\ that appear in the equations 

S(a H r,y/\ I) ,~ 

V ~~ * J, — ( ' ) 

a #oV I 1 - - | 

where S(x) = sinhx in an open universe and S(x) = sin a; in a closed universe. 



B. The Ostriker-Vishniac effect 



The reionization of the Universe is one of the most important physical processes that took place in the early 
universe (see |2fil2l| ). The most accepted sources for reionization, which requires a source of ultra-violet photons, are 
an early generation of massive stars formed in dwarf galaxies or an early generation of quasars/ AGNs in galaxies. In the 
currently favored adiabatic CDM class of models for structure formation, reionization is expected to occur in the range 
8 < z r < 30. Measurements of the CMB anisotropies on sub-degree scales |2^,[2^,[24| have been used to put an upper 
bound on the reionization redsfhit of z r ~ 30 |25l l26| . Very recently, u sing polarization and temperature anisotropies 
of the CMB, WMAP has placed a fairly model-independent constrain [23 on the optical depth to electron scattering 
of r = 0.17 ± 0.04 at 68% CL which translates into z r = 17 ± 3 for instant reionizatio n I2II 12^1 . Interestingly, the 
measurement of an increase of the neutral fractions with redshift in high-z quasar spectra |3Ctl3l| and a first detection 
of the 'Gunn-Peterson trough' [33] in a quasar spectra at z = 6.28 by SLOAN .33( point to a reionization redshift 
of z r ~ 6 34], in disagreement with WMAP results. However, even a fraction of neutral hydrogen as small as 0.1% 
in the IGM could explain the result due to the large cross-section to Ly-a photons. Together with the results from 
WMAP, high-z quasar measurements indicate that the reionization history is more complex than previously thought 
and attempts are being made to fully understand it [35L 1361 133 • 

Though extensive analytical (for a complete derivation see |l2L[l3| ) and numerical studies (see references below) have 
been done to try to quantitatively understand many of the effects originated by reionization on the CMB, accuracy 
is difficult to reach and uncertainties still remain. Reionization will leave multiple distinctive imprints on the CMB 
anisotropies by bringing the CMB photons and the free moving electrons into scattering contact again. Through that 
window, the low-z period of the universe evolution can be probed experimentally in more detail from the appearance of 
the first sources of ionization to the formation of the observable present large-scale structure. Studies have been done 
on the calculation of the contributions to the power spectrum of the CMB due to ionization induced effects like the 
Doppler effect on large angular scales I3SII3II. the thermal SZ effect and its kinetic analog l4pL l4lL 143 . | 4^. 14^ . I45L I46L |4^tT| . 
the Inhomogeneous Reionization |48l 14911501 IBH 152^ and the OV effect on smaller scales |l4l l53l l54| . Enlightening 
comparative studies between different effects can also be found |55| . 

As ionization effects introduce non-Gaussianities in the anisotro pies. further studies were done on the calculation of 
their possible contributions to the bispectrum. Many authors [56|, |53, IH, |53, llHl investigated contributions to mixed 
bispectra due to couplings between lensing effects, Integrated Sachs- Wolfe effect, thermal SZ and Doppler effects, 
such as the OV effect. The trispectrum of ionization secondaries |5!| hasn't been explored very much. No one has 
addressed the calculation of the pure bispectrum and trispectrum of the OV effect until now. 

In the linear regime, for the power spectrum, the dominant small-angular scale contribution from reionization was 
found to be the OV effect hj, Ml. It arises from the second order modulation of the Doppler effect by density 
fluctuations which affect the probability of scattering. Because of its density weighting, it peaks at small angular 
scales, typically arcminute scales in ACDM models, and should produce [iK anisotropies. Its contributions to the 
temperature fluctuations along the line of sight can be written in the manner of Jaffe and Kamionkowski (JK) |54| 

5f (0) = - £° drj g( V ) d.p(9rj, 77) (8) 

where p(#?7, 77) = v(0rj, r))5(9rj, 77) and g is the visibility function given by 

= a(y)aTnM e -r M = 0-138 O t ft + 2 ^ 
c c 

which gives the probability of scattering. The prefactor 0.138 is obtained assuming that all the baryons are in the 
form of protons (if we use the fact that the mass fraction of helium is 25% then one should multiply it by 7/8). 
The visibility function is normalized such that g{rf)drj = 1 — e~ Tr where r r is the optical depth to the surface 
of last scattering at recombination. Note that g is only dependent on time, and not on position for the OV effect. 
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The optical depth is given by t(t}) = L cdrj' or n e (?/). Also 5(9r],r]) and v(9r],T]) are the density contrast and bulk 
velocity along the line of sight, n e (v) is the mean electron density given by n e (r)) = p c x e (ri) (1 + z) 3 /m p , cr T is the 
Thomson cross-section, x e is the ionization fraction and m p the proton mass. We assume that the visibility function 
is approximated as a Gaussian in conformal time 

9(V)= /0 n - e^^^. (10) 

Following JK |54| , we choose a coordinate system such that 9 represents a three-dimensional unit vector along the 
line of sight and 6 refers to a two-dimensional unit vector in the plane perpendicular to the line of sight. So we will 
have = (9\, 6*2,0) and 9 = (61,62, y/l — 6\ — 9%) — (61,62, 1) where this approximation arises from the small-scale 
nature of the effect. Bold letters represent three-dimensional vectors. 

The OV is a small-angle effect so we can work under the fiat-sky approximation and expand the temperature 
perturbations in Fourier space 

^(K) = - r°dV9(v) I d 2 6 f-^-e.p^e^-^ (11) 



T v ' J IJK " J J (2k) 

where q = (q x , q y , q z ), K = (K x , K y , 0) and 

_ , , ia(n)i>D f d 3 k , . / q - k k \ 

p(q ' = I - k) (j^kF + w ) (12) 

is the Fourier transform of p(6rj, if) (see JK |H3|). D and D depend on r\. We made use of the continuity equation in 
Fourier space JSJ and of the linear evolution of the density field. 



C. Statistical properties of a general field 



The statistical properties of a field can be characterized by the n-point correlation functions in real space or by the 
n-point spectra in Fourier space. If the field is Gaussian in nature, like the primordial density fluctuations field in the 
current favored inflationary cosmology, the connected part of the n-point functions disappears for n > 2. The non-zero 
(even-n)-point correlation functions can be expressed with the 2-point correlation function. As a result, a Gaussian 
distribution is completely described by the two-point correlation function, or power spectrum, and any non-Gaussian 
field will be detectable by measuring the connected part of its n-point correlation function. 

If we consider a general statistically homogeneous and isotropic 2-dimensional field p with zero mean, its power 
spectrum P, bispectrum B and trispectrum T are defined by the following equations in the appropriate Fourier 
convention: 

<p(k 1 )p(k 2 )> = (2ir) 2 P(k 1 )S 2 D (k 1 + k 2 ) 
< p(ki)p\k 2 )p(k 3 ) > c = (2Tr) 2 B(k 1 ,k 2 ,k 3 )S 2 D (k 1 + k 2 + k 3 ) (13) 
<p(k 1 )p(k 2 )p(k 3 )p(k A )> c = (2TT) 2 T(ki,k 2 ,k 3 ,k A )5 2 D (k 1 +k 2 + k ? , + k±) 

where the subscript c stands for connected. The OV effect being a secondary effect will introduce non-Gaussianities 
in the original primordial Gaussian distributed temperature fluctuations. As a consequence, contributions to its 
bispectrum and trispectrum are expected. 



D. Dominant contributions among the statistics of the OV effect 



Combining the homogeneous theory of turbulence with the Limber approximation enables one to infer the domi- 
nant contribution among n-point statistics of an isotropic and homogeneous vector-like field effect whose statistical 
properties vary slowly in time. In particular, we can apply this to the OV effect. In short, the theory of homogeneous 
turbulence shows how to build invariant spectral tensors of arbitrary order, corresponding to expectation values of 
arbitrary products of statistically homogeneous vector fields. It is based on techniques proposed in the area of homo- 
geneous turbulence in the 1940s by Robertson [6l]|. Relying on this theory, all expectation values of an odd product of 
an isotropic 3-dimensional vector field p(q) with q = (q x , q y , q z ) must be proportional to at least one of the q vectors, 
contrary to the expectation values of even products. 



5 



Because of the Limber approximation, extended to higher-statistics in Appendix A, which states that the only 
contributions to the projected correlation function on the sky come from the Fourier modes perpendicular to the line 
of sight of the angular correlation function, all the qi z terms tend to be suppressed. There will be different levels of 
suppression depending on the order of the qi z dependence of our statistics. 

Combining these two results, we can conclude that even correlation functions of the OV effect dominate over odd 
correlation functions making the trispectrum a much more sensitive statistics than the bispectrum. Also, we expect 
the correlation functions to obey the homogeneous and isotropy theory fully and thus to be the result of contributions 
of different orders in the qi Z terms. We have developed a method which permits to calculate the dominant contribution, 
under the Limber approximation. 

This is a characteristic of all effects physically described by an isotropic vector field and can thus be useful for 
other studies. As noted previously by Scannapieco |63 |. the alternation of dominant /subdominant/dominant higher 
order correlation functions provides a unique signal distinguishing the OV effect from other non- vector like secondary 
anisotropics and consequently can enable one to disentangle it from other contributions at similar angular scales. 



E. Signal-to-noise 



A fundamental issue is to know how well we can separate the OV signal, which is non-Gaussian, from the Gaussian 
signal, noise and foregrounds which are always present in a measurement from an experiment. A way of quantifying 
this detection is to calculate the x 2 statistics (as in [HilEiil)- To do so, one needs to calculate the Fisher information 
matrix Fy (for a good review see ^ we think of the data x as a random variable with a likelihood function 

T(x; 0) where 9 is a vector of model parameters, the Fisher information matrix is defined as 

*._/2«teffl\ (14) 



dBtdOj 

By a very powerful theorem, called the Cramer- Rao inequality, it was shown |65l lol| that the variance of any unbiased 
estimator of a certain parameter in a model, can not be less than (F^ 1 )^. As the signal calculated is expected to 
be rather small, we are interested in estimating its overall detectability as in [HfiL IH7I Therefore, we assume that 
the form of our model 9 (in our case the bispectrum and the trispectrum) is correct and that the only interesting 
parameter is its amplitude A, where the true value of A = 1. Then the Cramer-Rao inequality tells us that the 
variance of the measurement of A is no less than a 2 (A) = (F~ 1 )aa and we define the x 2 statistics as 

x 2 -(§) 2 = ^ = (^W- (is) 

The calculation of the Fisher matrix (F)aa of the statistics of the OV effect involves the calculation of the contribution 
of the noise to the power spectrum C" mse , as it will be shown. The noise depends on the experiment characteristics. 

We consider a hypothetical experiment which maps a fraction of the sky f s k y with a Gaussian beam with full width 
at half maximum 9 f w hm and pixel noise a p = s/ y/t P i X , where s is the detector sensitivity and t P i X is the time spent 
observing each pixel. We use the inverse weight per solid angle, w = [a p 9 f w hm/To) 7 hi order to have a measure 
which is independent of the pixel-size [Hl^]. To = 2.73K is the CMB thermodynamic temperature. If only a fraction 
fsky of the sky is mapped, treating the pixel noise as Gaussian and ignoring any correlations between pixels a good 
estimate of the C£ oise H HI 111 is 

<™ e = hkyw- 1 e a > i{t+ V (16) 

where a, in radians, is the width of the beam if we assume it has a Gaussian profile. It is related to 9f w hmi m 
arcminutes, by — \/8ln2 9f w h m x 7r/10800. Note that if an experiment maps the full sky and then a fraction 
1 — fsky is subtracted, one should not multiply w^ 1 by f s k y (case of MAP and Planck). Hence we can estimate C™ mse 
for any experiment with characteristic f s ky, 9f w h m , s and t P i X . 

For the precise cases of MAP (renamed WMAP recently) and Planck, for which we used the specifications in tabled 
we need to take into account their multi-frequency coverage with different characteristics. The C™ mse is then defined 
as 57] 

(17) 



where the sum runs over all channels of the experiment and v is the frequency of the channel. 
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MAP 


Planck 


v (GHz) 


41 


61 95 


100 143 217 


Ofwhm (arcm) 


31.8 


21.0 13.8 


10.7 8.0 5.5 


Op OK) 


19.8 


30.0 45.6 


4.6 5.4 11.7 


/sky 


0.80 


0.80 



TABLE I: Experimental parameters for (W)MAP and Planck. 

F. Power spectrum of the OV effect 

1. Linear power spectrum 

In Fourier space, the flat-sky power spectrum of the OV effect is related (Eq. (|13|) ) to the following two-point 
expectation value of the OV temperature field perturbation AT/T 

AT - AT - 1 /" )0 f Vo f f 

< -^{Kij-^iKi) > = - I d mg ( m ) I d mg ( m ) I d 2 e, I d 2 e 2 



10 Jo 
d 3 qi f d 3 q 2 A 



(2tt) 3 J (2tt) 3 



71 i 



hj ■(< Pi(qi,m)Pj(Q.2,m) > + < Pi(<to,m)pj(<ix,'m) >) 



e i(Ki.6i— rn_q-i_§i) e i(K 2 .82— V212-82) (18) 

where p is defined as in Eq. 1)12(1 . Many authors have derived the expression for this statistics dUlIilElEl we 
use the JK formalism but a different technique which will be useful in what follows. 

As we see, this expression involves a double integration in time, angle and wavenumber, being numerically long to 
evaluate. It is useful to note that as the statistical properties of the field p vary slowly in time and as the OV effect is 
a small-angle effect we can employ the Limber approximation (see Al) to considerably simplify our derivations. We 
stress here that we are allowed to use this approximation as the OV effect takes place at sufficiently high I, where the 
difference between the approximation and the integral is very small. 

As the two permutations < pipj > are symmetric due to statistical homogeneity, we only consider the first one and 
multiply the result by 2. Using Eq. (|12fl for p and the Wick theorem for the Gaussian 3-dimensional density field 
correlation function which states 

< 5(ki)5( qi - k 1 )o(k 2 )o(q 2 - k 2 ) > = 

(2^) 6 F(fc 1 )P(| qi - k x IX^kx + q 2 - k 2 )<%( qi + k 2 - ki) + $(ki + k 2 )6 3 D (qi - k x + qa - k 2 )) (19) 

where P(k) is the power spectrum of density perturbations, we obtain two non-zero terms for < piPj > which can be 
written as 

<Pi(qi,J7ife(q2,r? 2 ) >= -2 GfaJGfo) i^(qi) ^(qi + qa) (20) 
with the time dependence functions gathered in G(rf) = I %a J^S ) and the general tensorial functions F a p given by 



= J d 3 K'P{a)P{b) + + (21) 

where a = K' and b = q^ — K'. We could now replace the two previous expressions directly in Eq. (|18|) but, for 
clarity purposes as it will become obvious soon, we refrain from doing so and instead we keep working with F a p. 

Indeed, in the small-sky approximation, for which the unit vector 9 ~ (0, 0, 1), we can contract the 6's of expression 
(|18H with the vectors a and b of the last expression (|21|l such that we are left with the line of sight components of a 
and b. We can thus define a new scalar function F such that 

F(q,) = 9 la § 2p F aP * J d 3 K'P{a)P(b) + = J d 3 K'P(a)P(b) (jf + ^ + |) . (22) 

The interesting step that follows is to expand this function in a z = K' z and b z = g JZ — K' z 

F („,, . / W W ) (*? (A + 1 - ^) + K - ^) + (4 )) . (23) 
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Expressions of this type will occur in the following and they illustrate the previous discussion in Sec. Ill Dl By looking 
at the integral, we see that the K' z integration for odd products of K' z is zero as the dependence on K' z of the terms 
with a and b is even. Therefore we are left with non-zero contributions from the first term (with K' z 2 ) and from the 
last term (with K' z ). By applying the Limber approximation we can also infer that the dominant contribution has to 
come from the K' z term as it has no explicit dependency on qi Z . As we know, in the Limber approximation, Fourier 
modes parallel to the line of sight (terms on g, z ) tend to be suppressed. We are then left with 2 terms which we 
will label dominant and sub-dominant terms depending on the order of their cancellation. We are interested in the 
dominant one. This cancellation was forecasted in Sec. Ill Dl in view of the homogeneous turbulence theory. Indeed, 
the power spectrum of the OV effect is expected to have contributions only from terms with no dependency on qi z 
and with a qi Z qj Z type dependency. 

After some straightforward algebra, in the Limber approximation framework, we find for the dominant contribution 

F(*)-/^P(«)P(5)(«?(^ + ^-^)) 

= -2nqi / dyi / dfiP(q i yi)P(q i y 2 ) 3 (24) 

Jo J-i Vi 

where we have performed a spherical coordinate transform such that fi — qi.K 1 , a = y\qi and b = y 2 qi with y 2 — 
\Jl + y\ — 2yifi. To obtain the components of K' z in the chosen coordinate system, we used the Limber approximation 
to assume qi Z ~ 0, such that K' z = K'y/l — /i 2 . This assumption preserves our dominant term but suppresses any sub- 
dominant term that could naturally arise when calculating the integral. As a consequence, in Eq. \2'.i\i. the dominant 
term (in q® z ) may contain hidden contributions to the subdominant term (in qf z ). That this indeed is the case can be 
understood by a very simple reasoning. Consider Eq. I|22l) . It is easy to show that the terms in a 2 and b 2 z give identical 
contributions to the integral, such that if we calculate twice the integral in a 2 , we should obtain the same result at 
the end, i.e., various terms depending on different orders in qi Z . By doing this, our term in g 2 z present in Eq. iL'-il) 
simply disappears. We might then expect it to show up in the integral calculation. But, most interestingly, when 
performing the calculations as previously using the Limber approximation, we obtain the same result, i.e., Eq. I|24|l . 
As a consequence, by imposing qi z ~ in this example, we are in fact suppressing a sub-dominant term we know 
should be present. In conclusion, there are more terms contributing to the sub-dominant power spectrum than the 
one present in Eq. 112. "ill and these aren't so easily calculated. Hence we neglect all the sub-dominant terms in the 
Limber approximation. In the following, for the bispectrum and the trispectrum calculations, similar problems are 
present but they will not affect the lowest order terms in qi z and even order in K' zl which are the dominant terms of 
interest to us. 

We can finally combine Eqs. (|18|l . I|20|) and 124|) . Applying Kaiser's method and the Limber approximation as 
described in the Appendix, we obtain the well-known expression for the dominant contribution to the linear OV 
power spectrum 

P ?ol(K) = ^ jf ^nv)} 2 (j^J S(K/ V )d V (25) 

where 

S(k) =kJ o °°d yi jf * dnP(ky 1 )P(k/l + y\- 2y^) { \~ fy^^y ■ ( 26 ) 

The Limber approximation reduced the dimension of the integral from 6 to 3 and helped to find an easier numerical 
and analytical expression. We note here that we obtain a difference of a factor 1/2 compared to JK, a discrepancy 
pointed out by them when comparing to previous work. 

For illustration we show plots of the linear dominant contribution of the OV power spectrum in figure Q for the 
fiducial ACDM model assuming z r = 8 and z r — 17. The correspondence between full-sky multipole moments Ct 
and the flat-sky Fourier space P(K) is straightforward Ci — P(K — £). As expected, the z r = 17 scenario (keeping 
x e = 1) increases the amplitude of the power spectrum, due to the rise of the optical depth, and shifts its peak 
towards smaller angular scales. Numerically, we find for the amplitude of the power spectrum the approximate scaling 
dependence C«~5oo — 7.5 x 10~ 18 x 2 log ' 4 (1 + z r ). We advise the reader to consult JK paper |5J| for the impact of 
changing the cosmological parameters (or the reionization history) on the power and peak of the effect as well as for 
the important detectability issues. The cosmological dependence applies as well for the higher-order statistics. 
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2. Nonlinear extension: the kSZ effect power spectrum 

This extension was calculated previously |69| and we present it for the dominant power spectrum for the sake of 
consistency. The kinetic Sunyaev-Zel'dovich (kSZ) effect results from a Doppler effect suffered by the CMB pho tons as 
they travel through large scale structures emerged in a bulk flow. As pointed out by Hu [6!j among others |4ll l5q . in 
adiabatic CDM cosmologies, nonlinearities only affect the density field below the coherence scale of the bulk velocities 
and so the nonlinear density held is uncorrelated with the large-scale velocity field, which remains linear. Hence, the 
mild nonlinear extension of the density contribution in the OV effect expressions naturally becomes the kSZ effect 
from large scale structures. 

As the result of this, we can use the previous calculations of the OV effect power spectrum to obtain a similar 
expression for the kSZ by introducing a nonlinear correction in the density field contribution. Following Hu [6flj | 
approach, we replace the linear density power spectrum with its nonlinear analogue but leave the contribution from 
the velocity power spectrum the same 

where P NL stands for the nonlinear CDM power spectrum and where the mode coupling integral S was calculated 
previously (Eq. I|26|)) under linear theory. This expression includes both the OV and the kSZ effects. 

To calculate the nonlinear power spectrum, we assume that the baryonic gas traces the dark matter |69| . We follow 
Hamilton et. al. [7(j who presented a scaling relation for the correlation function in the nonlinear regime that was 
generalized to its Fourier analog by Peacock and Dodds The basic hypothesis is that nonlinear fluctuations on 
a scale k arise from linear fluctuations on a larger scale 

fc Un = [l + A| o (fe)]- 1 / 3 fc (28) 

so that there is a function relating the nonlinear and linear power spectra at these two scales 

A^(A ; ) = / NL [A[ (Un) (fc li „)] (29) 

which can be fit to simulations. We use the ff^ proposed expression for /nl- For stronger nonlinearities other 
corrections are necessary. See for example [T^- One also needs the relation A|(fc) = fc 3 /(27r 2 )P(fc). 

As discussed in W.Hu this estimative should be seen as an upper limit of the kSZ because on very small 
scales the gas pressure, unaccounted for, smooths the gas density as compared to the the dark matter density. The 
assumption that the baryonic gas traces the dark matter was shown to break down at multipoles £ > fO 4 [69l l74|. 
We show a plot of the power spectrum corrected for mild nonlinearities in figure for illustration of the nonlinear 
enhancement of the power spectrum at small scales for z r = 8 and z r = f 7. 

III. THE BISPECTRUM OF THE OV EFFECT 

1. Linear bispectrum 

In analogy with the power spectrum, the flat-sky bispectrum of the OV effect is connected by Eq. (|f 3|) to the 
following expectation value 

AT AT - AT - — 1 f Va f Vo f Vo f f f 

< -jriK^—iK^iKz) > = — J dmgim) J q dmgim) j dmgim) J d 2 e 1 J d 2 e 2 J d 2 e 3 

'd 3 q x f d 3 q 2 f d 3 q 3 ~ ~ ~ . 

hiVijVM (< Pi(qi,m)Pj(q2,f72)Mq3, W > +perm.) 



(2tt) 3 J (2tt) 3 J (2tt) 3 

e i(-Ki.0i — r/iqi.fli) e i(K 2 .62-r) 2 <l2-62) e i(K 3 .8 3 -r) 3 q 3 .0 3 ) 



(30) 



where the five permutations are with respect to the ordering (qi,q2,q3)- Proceeding in a similar way we did for the 
derivation of the power spectrum, this numerically heavy expression to integrate can be considerably simplified by 
using a generalization of the Limber equation to higher-order statistics (see Appendix). 

We start by calculating a simplified expression for the first permutation < Pipjpi >. At the end we generalize the 
results to include the total six permutations. Using expression (|f 2|l for p and the Wick theorem for the Gaussian 
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FIG. 1: The linear (label L) and nonlinear (label NL) OV power spectrum for the fiducial ACDM model. The dot dash lines 
correspond to z r = 8 and the solid lines to z r — 17. In both cases we assume Az r = 0.1(1 + z r ) . 



3-dimensional density field 6-point correlation function < <5(ki)£(qi — ki)J(k2)<5(q2 — k 2 )5(k 3 )(5(q3 — k 3 ) > we obtain 
Cf .Cf/3! = 15 terms for < PipjPi > of which 8 are non-zero. After some simple calculations, these 8 terms can be 
condensed into 

< Pi(<li,Vi)P](<l2,m)Pi(<l3,V3) >= 4G(77i)G(r/ 2 )G(77 3 ) [I^-j(qi, q2) + Fiji(<li, q.3)] 8n(<li + 02 + q3) (31) 
with the time dependence function G(rf) = f *f§^ J and the general tensorial functions F a p 1 given by 

W%. *) - J d*K>P(a)P(b)P(c) - ^) (| - %) - 5) (32) 

where a = K', b = K' — and c = K' + q 3 . We concentrate on F and combine all the terms in Eq. (|31|l when 
appropriate. 

Again using the small-angle approximation for which 9 ~ (0,0, 1), we can contract the 0's with the vectors a, b 
and c such that we are left with the line of sight components of a, b and c. We can thus define a new scalar function 
F such that 

F(q i)q ,) = hJ^F^ * / d*K'P(a)P(b)P(c) - |) (| - ( J - J) . (33) 
Expanding this expression in different orders in K' z will give us the different levels of contributions under the Limber 
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approximation to the bispectrum 

Ffa.q,) = / d 3 K'P(a)P(b)P(c) 



K * ( ** ( a 2 b A + a i b 2 b 2 c 4 + b 4 c 2 ) + H* ^ fl 2 c 4 + ^4 + fl 4 c 2 6 4 c 2 ) ) H 

2/1 1 \ 2 / 1 1 \ / 2 2 



111111 



z 2 fe 4 a 4 fe 2 a 2 c 4 6 2 c 4 a 4 c 2 6 4 c 2 
2 1 1 



K ' V lz [ a 2 6 4 6 4 c 2 J + ^ { 6 2 c 4 a 2 c 4 J + 9iz<ljz { b*c 2 c% 2 



-133 "S3 1- ( 34 ) 



& 4 c 2 ) + ^ 2 ( 6 2 t 

Again, and as detailed before in Sec. Ill Fl the integral in K' z for odd products of K' z is zero, and we are left with 
non-zero contributions from the second term (with K' z 2 ) and the last term (with K' z °). The dominant contribution will 
come from the K' 2 though it has a dependency on q z which, by the Limber approximation, tends to be suppressed 
compared to terms which don't have such dependency. This result again confirms the discussion in Sec. Ill Dl which 
predicts that the bispectrum of the OV effect has contributions only from terms with a dependency on qi z or on 
products of 3 q z . 

We are interested on the dominant K' 2 term which can be written as 

F(qi,qj) = gw/i(q»,qi) +&*/2(q<,qj) (35) 

where 

- / d^K'P{a)P{b)P(c)K' 2 + - ^ + (36) 

and 

/a(*,qi) = / d 3 K'P(a)P(b)P(c)K' 2 (-^ + ^ + -L. . ( 37 ) 

If we were to replace directly these last three expressions in l|30|l and then apply the Limber approximation we 
would have a direct cancellation. We wish to obtain an expression non-null which enhances the dominant character, as 
compared to the other sub-dominant terms, of this second-order contribution. To do so, we use again the small-angle 
approximation and consider 9i ~ (0, 0, 1) to write 

F(q 2 ,qj) = 9i.qifi(qi,qj)+§j.q j f 2 (qi,qj). (38) 

Now we replace this particular equation in 131fl and then in (|3UH and finally integrate once by parts the time de- 
pendence. This is the step that allows us to keep this term, which was expected to be suppressed under the Limber 
approximation. Then we follow Kaiser's method and apply the Limber approximation, and using (|13fl we obtain the 
expression for one of the six possible permutation terms of the bispectrum 

i r Va h 2 r)h 

B° v m (K u K 2 ,Ks) = t^3 7 o d V U^^) + (39) 

where h(r}) = (s^^p) an d = /i + h(&,qj)- 

In order to simplify the functions f± and fy, we need to find an explicit relation between a, b, c and q^, q_, and 
K'. To do so, we express both K' and a, b, c in the basis (e z , q~i, <jj) where we remind the reader that q a ~ (q a ±, 0) 
under the Limber approximation. Again, as for the power spectrum (see discussion in Sec. Ill F^ . this suppresses any 
sub-dominant contribution that could arise in the process but leaves our dominant contribution intact. We obtain 

/(qi,q?) = /(I q* I, I q? 

= fdKfduJ d V VT-^P(a)P(b)P( c)K 2 (-!_ + -I_ - - + _L + _L) (40 ) 
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where a = (K' z 2 +u 2 + v 2 + 2uvfi) 1 ^ 2 , b = (K' 2 + (u-q i ) 2 +v 2 + 2(u-q l )vfi) 1 ^ 2 , c = (K' z 2 + (v+q j f+u 2 +2(v+q j )un) 1/2 
and /i = cji.qj. As this expression corresponds to only one of the permutation terms of expression (|3Ufl . we need to 
include the other 5 to obtain the final expression for the dominant flat-sky bispectrum 

B$Z l (KuK 2 ,K 3 ) = -t 1 £ f(^A) (41) 

8^ 3 Jo t dr, . . \ r, n ) 

where / is given by (|40J) . We thus reduced our initial expression with twelve integrations to a four-dimensional 
integral, which can be numerically calculated for a chosen configuration of the wavenumbers K . 

Note that we obtain the first order time derivative of h. It involves one single derivation, which indicates the <& z 
cancellation. As h is smooth and slowly varying, the contribution from this term should be very small. Hence, the 
Limber cancellation at small scales reveals itself in the derivatives of the time dependent functions. 



2. Nonlinear extension: the kSZ effect bispectrum 



We follow the same approach as for the power spectrum in Sec. Ill Fl and apply it to the dominant contribution. As 
we have three CDM power spectrum showing up in the expression for the bispectrum, which we need to divide among 
linear/nonlinear contributions from density/ velocity contributions we will assume the following bispectrum effect 




3. The Signal-to-noise 

To define the \ 2 statistics of Eq. I)15|l. we need to calculate the likelihood L of observing the bispectrum elements 
Bp = Bg^^ given the true parameters p, and calculate the Fisher matrix as defined in Eq. (|14fl 

Assuming that the likelihood is Gaussian, we follow Cooray and Hu approach to calculate the \ 2 statistics 



x 2 



where B^ 1 i 2 i 3 is the angular averaged bispectrum defined on the sphere, f s k v represents the reduction in signal-to-noise 
due to incomplete sky coverage and vf t t comes from the covariance matrix of the angular averaged bispectrum 
assuming a nearly Gaussian bispectrum and full-sky coverage jH^, F(| 



°W S = C C C Z C Z I 1 + Wi + + 5 6 D & + 4) + 52, (* 3 + h) + 2SU*i + + 4)] (45) 

where C| ot stands for the sum of the power spectra of the primary cosmic signal, the thermal SZ (thSZ) effect 
which contributes significantly at the scales of interest, the linear OV effect, the detector noise and the foregrounds 
respectively 

Qtot _j_ QthSZ _j_ qOV~L _|_ Q-noise _|_ q foregrounds (46) 

The thermal SZ effect was taken from [7^ and was calculated semi-analytically at 30 GHZ for a normalization factor 
(Ts = 0.9. We include the linear OV effect contribution (see Sec. Ill F|l although its amplitude is small as compared 
with the primary and the thermal SZ signals at the scales considered in this work. The primary cosmic signal was 
computed with CMBFAST [73 and we will not consider any foreground. Indeed, for the case of MAP and Planck, 
studies indicate that the total Ci should increase by 10% maximum However, caution is required as this result 
assumes that foregrounds have a Gaussian distribution. The foregrounds contribution to the higher-order statistics 
could in fact be our main obstacle in measuring non-Gaussian effects and very little is known about them. For the 
noise power spectrum we use Eqs. (|16H17f) applied to MAP and Planck. 
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We will study the contributions to the \ 2 P er log interval in I, It gives us more sensibility on the angular scales of 
stronger detectability, depending of the effects considered. This also enables us to directly compare our results with 
Cooray and Hu [H3| . We calculate the total S/N for each experiment by integrating over the multipole £. 

A useful relation between the flat-sky bispectrum and the spherical harmonic angular averaged bispectrum [(| is 



/ (2l 1 + l)(2l 2 + l)(2l 3 + l) (h £2 M H ,„ . K . K f s U7) 
Vi&h = V |^ I j ' 1 = u 2 = ( ' 2 ' J<3 = ( 47 ) 

Since the Wigner-3j vanishes if l\ + £2 + £3 — odd, the full-sky bispectrum can only be estimated for even terms. 



4- Results 



We are concerned with the overall detectability of the dominant term, previously calculated. Hence we choose the 
simplest of the possible configurations, highly localized in Fourier space, K\ — K2 = K3 = I, for which the flat-sky 
dominant bispectrum becomes 



,<>\ ■ 3* r V0 , h 2 dh 



B» dL{ K = £) = - s ] Q *, ¥JRj f[-) (48) 

where / is defined by Eq. I|40|) . The nonlinear analog follows from the previous equation. 

With this configuration, depending only on the multipole £, there is a simple way of calculating an estimate of 
the order of magnitude of the (S/N) 2 per bin of £ were we to include all the remaining £ modes f6j of the full-sky 
bispectrum. Indeed in the Eq. (|44|l . we see that the number of modes contributing to the (S/N) 2 per bin oi£ increases 

( £ £ £\ 2 

as £ 2 and in Eq. 147|) ^ 3 ( q / mcreases as 0-36 £ so 

dx 2 , p Bj f2f3 (£ £ £\ 2 B 2 (£) 3 B 2 (K = £) 

-jl~fskyt ^T~/.*y* * {000) —~0-^fsky£ -2 ■ (49) 

The a 2 is calculated using Eq. (|4*5|l for all the £s equal and B(£) using Eq. (|4*%|l . 

We show the plots for the linear and the nonlinear flat-sky OV bispectrum in figure J2J for z r = 8 and z r = 17 in 
the context of our fiducial cosmological model. For z r — 8 the peak of the effect occurs around multipole I ~ 500 with 
an amplitude of l 3 B(£)/(2ir) ~ 1.2 x 10~ 23 whereas for z r — 17 the peak takes place at a higher multipole of I ~ 700 
with a stronger amplitude of £ 3 B(£)/(2ir) ~ 3 x 10 -23 , as expected. We find numerically for the amplitude of the 
bispectrum the approximate scaling dependence with reionization history B(£ ~ 50) ~ 1.7 x 10 -29 x\ log(l + z r ). 
The dependence with z r is neither in agreement with the one obtained for the power spectrum (see Sec. Ill F|) nor 
with the one obtained for the trispectrum (see Sec. liVfl . This can be explained by looking at the bispectrum Eq. I|48|) 
where we have the first order time derivative of h, contrary to the corresponding expressions for the power spectrum 
(Eq. ) an d the trispectrum (Eq. (|65|0 . This derivative of h introduces a stronger scaling relation with z r for the 
amplitude of the bispectrum as compared to the one for the other two statistics. For both reionization scenarios, 
the most interesting feature is the rapid drop of power after the peak which is not observed for the power spectrum 
(see figure JIJ), which is a direct consequence of the Limber cancellation at small angular scales. The fact that the 
bispectrum is not considerable enhanced by nonlinearities, which take place at small angular scales, is also the result 
of the effect peaking at intermediate multipoles. The second, but expected, result is the low overall amplitude of the 
effect. The z r = 8 case (which we can easily compare with previous results in the literature) is lower by more than 10 
orders of magnitude than the lensing couplings presented by Cooray and Hu [53] , and by a few orders of magnitude 
than most of the OV couplings involving the SZ effect though it should be comparable to the OV couplings involving 
the Doppler or ISW/SW effects. Despite the fact that such OV couplings suffer the Limber cancellation as well (they 
correspond to expectation values of an odd product of a vector field), this cancellation can be counterbalanced by 
the coupling to higher amplitude effects, like the SZ effect, and by the matching in redshift between the density 
and velocity fields of the OV effect and the secondaries to which is couples. This is the case of the hybrid coupling 
ISW-SZ-OV presented by jH^, that has the largest signal of all secondaries that couple with OV. 

A more revealing quantity than the power of the effect is the signal-to- noise ratio. We plot in figure J5J), the 
estimated contributions to \ 2 per log interval in £ of the linear full-sky bispectrum had we included all the modes 
for MAP, Planck [l| and no instrumental noise. We show the results for z r = 8 and z r — 17. The experimental 
specifications can be found in table [I] We found no necessity of plotting the corresponding d\ 2 /d£ for the nonlinear 
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FIG. 2: Left panel — Linear (label L) flat-sky bispectrum of the OV effect and its nonlinear extension (label NL). Right panel 
— Contribution to \ 2 P er log interval in £ for the OV full-sky linear bispectrum with no instrumental noise (top), Planck 
noise (middle) and MAP noise (bottom) included in the variance. We used the specifications in the tableU] All the plots were 
calculated for the fiducial ACDM model. The dot dash lines correspond to z r = 8 and the solid lines to z r = 17. We assumed 
Az r — 0.1(1 + z r ). The total S/N for the OV full-sky linear bispectrum for MAP, Planck and a perfect experiment respectively 
are: 5.0 x 10~ 5 , 1.2 x 10~ 3 and 4.4 x 10~ 3 assuming z r = 8, and 3.7 x 10~ 5 , 2.0 x 10~ 3 and 1.4 x 10~ 2 assuming z r = 17. 



bispectra due to its similarity to the linear bispectra. The structure in the d\ 2 /d£ arises mainly from the structure 
of the CMB primary power spectrum at £ > 200. We point out that the continuous rise at I ~ 3000 for a perfect 
experiment is due to the fact that the signal decreases slower than the contributions to the Ci from primary linear 
OV and thermal SZ anisotropies up to these scales. Though it is not explicit in the figure, the d\ 2 /d£ is zero when 
31 is odd. As we can see, considering thermal SZ contributions to the noise evaluated at 30 GHz, the signal-to-noise 
of the OV bispectrum is very small, even for a perfect experiment, for which the total S/N ~ 4.4 x 10~ 3 for z r = 8 
and S/N ~ 1.4 x 10~ 2 for z r = 17. These values can be compared to the value S/N ~ 1.7 obtained by Cooray and 
Hu [57J for the bispectrum generated by the coupling ISW-SZ^-OV. Hence, even for early reionization, no detection 
is to be expected from future experiments unable to remove the thermal SZ effect. For the ideal case of a perfect 
multi-frequency experiment observing in the millimeter and sub-millimeter capable of subtracting all of the thermal 
SZ effect, the total signal-to-noise increases to S/N ~ 0.1 for z r = 17, which means that not even for this case are 
we to expect a detection. Our results for the S/N rely on various convenient assumptions required to considerably 
simplify the analytical expression for the estimate of the S/N (see sections Sec. Ill El and Sec. 1111 Sf> . This means 
that the values obtained for the S/N are most probably upper limits on the real values expected and should thus be 
considered accordingly. 
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IV. THE TRISPECTRUM OF THE OV EFFECT 



1. Linear trispectrum 

We follow the same procedure used for the power spectrum and the bispectrum. The trispectrum contribution to 
the temperature fluctuations is connected by Eq. i(T3|) to 

AT - AT - AT - A~T - 



dvigfa) / d mg ( m ) / dr) 3 g( m ) / d mg ( m ) / d z e l / d z e 2 / d% / d?e 



i 

24 Jo 



(27r)3 i(27r)3 JThrf J J^nf 9u92 i e ^ eim (< P<(<H« ^i)Pj(Q2) m)Pi{^i, m)Prn{<U,r)4) > + perm.) 

e l(/f 1 .0 1 -?) 1 q 1 .e i ) e l(X2-S2-J)2q2-62) e i(^3-S3-J73q3.e3) e i(if4.e4-774q4.e4) /gg\ 

The 24 total permutations arise from symmetries under permutation invariance. Again we concentrate on the first of 
the permutations and then generalize at the end. We obtain this time for the first permutation term of the previous 
Eq. < pipjpipm > by the Wick theorem applied to the Gaussian 3-dimensional density field 8-point correlation 
function Cf .C|.C|/4! = 105 terms of which 12 non-zero terms for the Gaussian contribution and 48 non-zero terms 
for the connected part which interests us here. Performing some exhaustive and systematic calculations these 48 
terms can be written in the following condensed form 

< Pi((ll,Vl)Pj((l2,V2)Pl(<l3,1l3)Prn((li,Vi) > = 4 G^G^G^G^) 

[Fimij (qi,q2,q3,q4) + ifym (qi,q4,q2,q3) + -Fkim(q3,q4,qi,q2) + 

Fijmi (q3, qi, q2. q4) + F jmH (q2, q 3 , qi, q4) + F jilm (q2, q 4 , qi, q 3 )] 
<5 3 (qi +q2 +qs +q4). (51) 

where G is defined in the previous section and the general tensorial function F a p 7 s is defined as 

F a @ 7 8(<li,<lj,<ll,<lm) = fa/3<ys(<li,<ij,<ll) + fa"//36(qi,<lj,<hn) (52) 

where a general f xy zw is given by 



(53) 



where a = K', b = K' — q.j, c = K' + q^ and d = K' + q^ + q;. Again we will work with a single / and at the end 
generalize the result. 

Applying the small-angle approximation again we can remove the tensorial dependence of our expression and 
contract the 9i with our vectors. We obtain the scalar / 

/(qi>q3,q;) = flia^^^/o^sCq^qj.q;) 

« / m)P{mc)m (I + I) (I + *) - *) (54, 

and the scalar F 

^(q»,qj,qj,qm) = /(qi,q?,qO + /(q 4 ,qj,q™)- (55) 

Expanding / in order of iC£ and keeping the non- nulls terms (see previous section), i.e. the terms which are even in 
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K' we are left with 



/(%,q,-,qj) = J d 3 KP(a)P(b)P(c)P(d) 



K' A 



1 
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4 2 2 1 2 1 2 1 



-QizQlz 



b 4 c 4 a 2 b 2 c 4 a 2 b 4 c 2 a 4 b 2 c 2 a 2 c?c 4 b 2 c 2 d 4 a 2 c 4 d 2 a 4 c 2 d 2 
2 2 2 1 1 2 2 



i d 4 a 2 c 2 d 4 a 2 b 4 d 2 a 4 b 2 d 2 b 2 c 4 d 2 b 4 c 2 d 2 a 2 b 2 c 2 d 2 

2 112 



I a 2 c 2 d 4 6 2 c 2 d 4 + fl 2 c 4 d 2 



6 2 c 4 d 2 a 4 c 2 d 2 6 4 c 2 d 2 a 2 & 2 c 2 d 2 



.3/1 1 \ , 2 2 I 1 1 



( fer?.?* 19 2,J2 (,4„2J2 + 



6 2 c 2 d 2 6 4 c 2 d 2 J \b 4 c 2 b 2 c 4 d 2 

(^h - v&p) + + (p^) } ] ■ (56) 

We concentrate on the dominant term which can be simplified using the same method applied to the calculation of 
the dominant term of the power spectrum. It can be written as 

= / *K>P(a)PQ>)P{c)PWK« (i - i) (i - i) (1 - I) (i - 1) . (57) 

To find an explicit relation between a, b, c, d and q^, q^, q; and K' we express both K' and a, b, c, d in the basis 
(e z , qt, qj). We should say here that using the Limber approximation will allow only two of our initial vectors among 
the set of four (qi,q2,q3,q4) to be independent. Indeed all parallel components to the line of sight will be negligible 
and the four vectors will be inside the same plane, the one perpendicular to the line of sight. This justifies the use 
of the basis chosen. We stress one time more that this eliminates any sub-dominant term that could show up (see 
discussion in Sec. Ill F|l . So in that basis our function / can be expressed in the following simplified way 



/(q^q^q;) = f(\ q 4 1, 1 q? 1,1 q* \Ai4jAi4iAj4i) 

i :i \ ( :i i \ ( i i \ ( l l 

c 2 a 



= IdK'Jdul dvVT^P(a)P(b)P(c)P(d)K? (--- --- -_-)(---( 



where a — q t -q^i P = Qi-Ql an d 7 = Qj-Ql' Also a — (K' 2 +u 2 + v 2 + 2uva) 1 / 2 , b — (K' 2 + (u — q t ) 2 +v 2 + 2(u — q^va) 1 / 2 , 
c = (K 12 + {v + q 3 ) 2 +u 2 + 2(v + q j )ua) 1 / 2 , d = (K 12 + {v + q 3 +y) 2 + (u + x) 2 + 2{v + q J +y)(u + x)a) 1 / 2 with y = qi-^E^j 
and x — qi ME^] ■ Our / on ly depends on the norms of its arguments and on the angles between their directions. We 
can combine these expressions to simplify the expression 1)5()[). Following Kaiser's method and proceeding with the 
Limber approximation, we obtain the expression for one of the 24 possible permutation terms of the trispectrum 

T dom {K u K,,K,,K,) = —sj Q dnl^j - r 

[FiKx/v, K 2 /r), K 3 , /r)Kih) + FiK./r,, K A /r), K 2 /v, K 3 /rf) + F(K 3 /r,, K 4 / V , K^r,, K 2 / v ) + 

F(K 3 /r], Ki/ V , K 2 / v , K^h) + F(K 2 /v, K 3 /v, Ki/rj, Ki/rj) + F(K 2 /v, K^h, Ki/v, K 3 /v) ] (59) 

where F is defined by Eq. I|55|l and the / by Eq. (|58|l . Finally, by including all the permutation terms in Eq. I|5U|I we 
obtain 



8tt3 J \ D 2 j tf \ V V V 



(60) 
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where / is given by 158|l . The power of the Limber approximation was to reduce an almost impossible integration to 
a 4-dimensional integral, which can be numerically calculated for a chosen configuration of the wavenumbers K. 

2. Nonlinear extension: the kSZ effect trispectrum 

We follow the same approach as for the power spectrum and bispectrum and again we calculate the nonlinear 
extension for the dominant contribution. As we have four CDM power spectrum showing up in the expression for 
the trispectrum, which we need to divide among linear/nonlinear contributions from density/ velocity contributions 
we will assume the following trispectrum statistic 

"•---'-^r^^)^,X,( J ») ,/ (T44)- <«> 



3. Signal-to-noise 

As pointed out by Zaldarriaga [j^ and later by Hu [6!| , the maximal signal-to-noise can be proven to be 

x 



N ) J sky 2^ 2^ 2L + 1 C tot C tot C tot C tot ^ 

' L £i<£ 2 <f3<^4 £l * 2 ^ f ' 4 



where T£ 1 ^ 2 ^ 3 ^ 4 (L) is one of the possible configurations of the full-sky trispectrum as defined in what follows. The 
covariance matrix used to obtain the Fisher matrix is calculated assuming full-sky coverage and l\ < £2 < £3 < £4 
by Komatsu J6(J in the weakly non-Gaussian limit. If l\ < £2 < £3 < £4 is not respected, the covariance would 
be distributed across many Ls and can lead to overestimates of the signal-to-noise By not respecting this last 
constrain in Eq. (|62(l we are calculating an upper limit of the S/N estimate. 

Concerning the equivalence between the full-sky and flat-sky formalisms, we follow Hu's Appendix [|| where it is 
argued that we can find a relation between the two formalisms by breaking up the trispectrum in the three possible 
combinations in each configuration defined by (£1, £2, £3, £4) 

2V a £ 3 i 4 = T(£ 1 £ 2 )(£ 3 £ 4 )(Li 2 ) + T(^ 3 )(£ 2 £ 4 )(Li3) + T(£ 1 £ 4 )(£ 3 £ 2 )(Li 4 ) (63) 

where Ly correspond to the side of the triangle of sides ij, Ly). Each of the T^ i i j ^g k g m )(Lij) is then related to 
the flat-sky equivalent by 

T&w^Lu) = ^±1^ + l){2£ 3 + 1)(24 + \){2£ m + 1) ( \ \ ^ ) ( ^ ^ ^ ) 

T((Ki = £ h Kj = £ j ), (K k = £ k , K m = £ m ){Uj)). (64) 

As for the bispectrum, the Wigner-3j vanishes if £i + £j + Lij = odd. So the full-sky trispectrum can only be estimated 
for even terms. 

4- Results 

We are interested in the numerical evaluation of the dominant contribution of the trispectrum as it may be of 
cosmological interest in the near future due to the next generation of experiments. We choose the trapezoidal 
configuration in Fourier space K\ = K2 = K3 = K4 = I with an angle 9 between two consecutive sides. This will give 
rise to the following flat-sky trispectrum 
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where e = cos(9) and / is defined by Eq. (|58|l . The various / correspond to the different specific configurations due 
to the ordering of the vectors Ks in the sum over the function / in Eq. (|60|) . Note that the analytical expression of 
Tdom 1S symmetric around 8 — tt/2 reflecting the symmetry of the configuration. The nonlinear analog follows from 
the previous equation. 

With this configuration which depends only on the multipole I and the angle 9 between two sides one can calculate 
an estimate of the order of magnitude of the (S/N) 2 (9) per bin of £ were we to include all the remaining t modes of 
the full-sky trispectrum. Indeed in the Eq. I(62|l . we see that the number of modes contributing to the (S/N) 2 per bin 
of £ increases as £ 3 so we will have the following relation for \ 2 m a bin of 1 

___!_) „, 3 f y I UU) | 2 

di Jsky ^ {2L l + l)(Cf t Y 1 ' 

where Tg(Li) is one of the 3 possible configurations of the full-sky trispectrum for a given multipole £ and angle 9 
which needs to be calculated from its flat-sky counterpart (Eq. (|65[l ) using Eqs. (|63[) and l|64f> . The C| ot is calculated 
as for the bispectrum (see Eq. ffify ). Again, f s k y represents the reduction in signal-to- noise due to incomplete sky 
coverage. 

We show the plots for the linear and the nonlinear OV flat-sky trispectrum for —0.95 < e < 0.00, assuming both 
z r = 8 and z r = 17, in figure J2J). As the trispectrum is symmetric around e = 0.00 we only plotted the negative es, 
but this choice was arbitrary. For z r = 8 the linear trispectrum peaks around multipole £ ~ 2.10 3 and has a maximum 
amplitude of £ 4 T '(£) / (2ir) ~ 2.4 x 10~ 27 , regardless of the configuration. Choosing an earlier reionization of z r = 17 
shifts the peak to £ ~ 3.10 3 and increases the overall amplitude of the trispectrum to l 4 T(£)/ (2ir) ~ 1.2 x 10 -26 . 
Numerically we find for the amplitude of the trispectrum the approximate scaling dependence with reionization history 
T{£ ~ 170) ~ 2.6 x 10~ 31 x\ log ' 8 (1 + z r ), in agreement with the relation found for the power spectrum. The effect 
of increasing the value of angle between two consecutive sides of the trapezoid is to increase the power at small scales 
(I > 10 4 ), such that only small scales are sensitive enough to probe different configurations. At those small scales, 
the power is maximum for e = 0.00 and then decreases as we decrease the angle. So at small angular scales the square 
configuration is the one contributing the most. We could have expected this behavior as the OV effect has a quite 
symmetric morphological signature, such that a more filamentary structure probed by the collapsed configuration of 
the trispectrum is not as likely. 

Here we do not observe the sharp drop in power at small scales due to the Limber cancellation which enables one 
to speak of the trispectrum of the kSZ effect. Indeed, concerning the nonlinear trispectrum, we observe an interesting 
feature. Contrary to the bispectrum and similarly to the power spectrum, the trispectrum is strongly affected by 
the weak nonlinear enhancement due to formation of structure at small scales. This enhancement has the power to 
broaden the shape, to shift the peak of the effect to smaller angular scales and to increase slightly the amplitude 
in a configuration dependent way. So (for z r — 8/z r — 17) we measure a higher and higher maximum amplitude 
(~ 2.4 — 7.2 x 10~ 27 /~ 1.9 — 4.8 x 10~ 26 ) at a smaller and smaller angle (£ ~ 1 — 4.10 4 ) when you go from the collapsed 
trapezoid to the square configuration. It is worth noticing that, whereas the amplitude of the nonlinear trispectra 
increases, the multipoles corresponding to the peak of the effect remain the same for both reionization scenarios. This 
is because the nonlinearities take place at the same instant in time for both reionization scenarios, leaving an imprint 
at the same characteristic angular scale. 

But the most important quantity is the signal-to-noise. We show in figure (@J) the estimated contributions to x 2 
per log interval in £ for the linear and the nonlinear OV full-sky trispectrum (z r = 8 and z r — 17) had we included 
all the modes for e = —0.95 with no instrumental noise, Planck noise and MAP noise included in the variance 
We chose e = —0.95, i.e. the collapsed configuration, because it generates the highest contributions to the full-sky 
trispectrum near the peak of the effect, contrary to the flat-sky trispectrum for which the square configuration was the 
one producing the highest amplitude. The square configuration continues to provide the strongest contribution at very 
small angular scales, but at larger angular scales the collapsed configuration dominates. This is due to the angular 
averaged factors relating flat- to full-sky trispectra (see Eq. 1)64(1 ). Though it is not explicit in the figure, the dx 2 /d\n£ 
is zero when 2£ + L is odd. Again, as for the bispectrum case, the structure in the d\ 2 1 d£ arises mainly from the 
structure of the CMB primary power spectrum at I > 200. The continuous rise at I ~ 3000 for a perfect experiment 
is due to the fact that the signal decreases slower than the contributions to the Cg from primary, OV and thermal 
SZ anisotropies up to these scales. Contrary to the common and most naive expectation, for Planck, an eventual 
detection is possible at arcminute scales. Indeed, the total S/N ~ 4.7 for z r = 17. Most of the contributions come 
from multipoles between £ = 1 — 2.10 3 . This probably is the most important conclusion of this work and illustrates our 
predictions. This result should be taken with caution as it corresponds to a very optimistic upper limit on the S/N 
(see sections Sec. Ill El and Sec. II V 3|) . Firstly, the use of the Fisher matrix formalism gives the minimum variance for 
our statistic (see section Sec. Ill Ell . Secondly, the x 2 method assumes that the form of the model is correct, which may 
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FIG. 3: Left panel — Linear and Right panel — nonlinear flat-sky trispectrum of the OV effect for geometrical configurations 
such that —0.95 < e < 0.00 in steps of 0.05. The amplitude of the trispectrum decreases as e decreases from 0.00 to —0.95. 
Because the power is symmetric in e around 0.00 we only plotted the negative es. All the plots were calculated for the fiducial 
ACDM model. The dot dash lines correspond to z T = 8 and the solid lines to z r = 17. We assumed Az r = 0.1(1 + z r ). 



not be the case. Thirdly, when calculating the covariance matrix of the Fisher matrix, two simplifying assumptions 
were used: that the main contribution to the covariance was Gaussian in nature and that we observed with full-sky 
coverage. Finally, other physical mechanisms, like for example lensing effects or the thermal SZ effect pi l59l l79l l80| . 
and unaccounted foregrounds [T^ can contribute to the trispectrum at this level and so the separability problem needs 
to be addressed in due time. Of course, uncertainties of roughly an order of magnitude in the modeling of the thermal 
SZ signal are also a source of error in our estimates. The forecasted ability of future multi-frequency experiments to 
remove most of the thermal SZ contributions would minimize these uncertainties and would much favor a detection. 
Last but not least, further progress in the implementation of optimal unbiased trispectrum estimators to probe such 
small scales and power is required. 



V. CONCLUSIONS 



Due to its strong predictive power, linear theory is a very sensitive probe of the early stages of the reionization 
history through the Ostriker and Vishniac effect. Analytical expressions for its correlation functions can be derived 
and their measurement would be of high value to our present knowledge of that still unclear epoch of the universe 
evolution. We have presented detailed calculations of the three Fourier statistics of interest of the OV effect, the power 
spectrum, the bispectrum and the trispectrum. For that purpose we have developed a new technique that allows one 
to obtain their dominant contributions under the Limber approximation framework. This method is applicable to the 
derivation of any statistics involving correlations among vector-like effects. It illustrates what was expected under 
statistical homogeneity and isotropy assumptions and the vector and small scale nature of the OV effect. We also 
evaluated numerically, as a function of scale and for a specific configuration (equilateral for the bispectrum and 
trapezoidal for the trispectrum), these statistics for a flat ACDM cosmology and two reionization scenarios. The first 
one is based on our pre-WMAP knowledge (z r — 8) and the second one takes into account the high values for the 
electron optical depth measured by WMAP (z r — 17). We numerically obtained approximate scaling relations for 
the amplitude of the OV statistics on the reionization history considered and found that the dependence is stronger 
on the ionization fraction than on the redshift of reionization. We have also studied their detectability in view of 
future satellite experiments. The alternation of dominant/subdominant /dominant higher order correlation functions 
was numerically shown. While the bispectrum is undetectable even by a perfect multi-frequency experiment capable 
of subtracting the thermal SZ contributions, the trispectrum could be measured by Planck or by interferometer 
experiments targeting arcminute scales with high sensitivity and for a sufficiently long period of time. This provides a 
unique signal distinguishing the OV effect from other non vector-like secondary anisotropics and could be useful when 
trying to separate different physical mechanisms imprinting themselves on these measurable statistics. One should 
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FIG. 4: Contribution to the x P er l°g interval in £ for e = —0.95 for the OV full-sky linear/nonlinear trispectrum with 
no instrumental noise (top), Planck noise (middle) and MAP noise (bottom) included in the variance. Again we used the 
specifications in the table H] All the plots were calculated for the fiducial ACDM model. The dot dash lines correspond to 
z r — 8 and the solid lines to z r — 17. We assumed Az r = 0.1(1 + z r ). The higher amplitudes for each of the experiments 
correspond to the contributions from the full-sky nonlinear trispectrum. The horizontal line at d\ 2 /dln(l) = 1 shows the 
minimum detection threshold. The total S/N for the OV full-sky linear trispectrum for MAP, Planck and a perfect experiment 
respectively are: 1.4 x 10 -3 , 1.1 and 27.4 assuming z r — 8, and 2.9 x 10 -3 , 4.7 and 149 assuming z r = 17. 



bear in mind that despite this useful characteristic signature, our results are quite optimistic, although encouraging, 
as they rely on various analytically helpful idealized assumptions, as described previously. Also, other contributions to 
the signal are to be expected at the arcminute level and thus further study of CMB small-scale secondary anisotropics 
and foregrounds contributions to the trispectrum is required and much justified. In order to obtain an upper limit 
on the possible kSZ contributions, we also extended our calculations to the mildly nonlinear regime. We found 
that, contrary to the bispectrum, there is a noticeable enhancement of the contributions of the trispectrum in a 
morphological dependent way and that this enhancement reflects itself on the calculations of the signal-to-noise. 
Hence nonlincarities are expected to enhance the even non-Gaussian signals produced by the OV effect and further 
complicate its disentanglement from inevitable model-dependent nonlinear effects arising from structure formation. 
Conversely, having a template of the OV effect can help in extracting the nonlinear contributions to reionization at 
those angular scales, providing another possible window on the complex physics associated with reionization. 
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APPENDIX A: LIMBER APPROXIMATION 

We review the fundamental steps of the limber approximation as used in the text. The Limber Eq. J^H describes 
the two-point statistics of a field which is the two-dimensional projection on the sky of a three-dimensional field whose 
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statistical properties vary slowly along the line of sight. The Fourier space analog of this result was calculated by 
Kaiser [82^ . Later it was further extended from the flat-sky approximation to an all-sky approach for spatially flat 
cosmologies by Hu and White [&| an d to open cosmologies by Hu Buchalter et. al. [13 derived the Fourier 

space analog of the Limber's equation for the bispectrum. 

Here we review the bispectrum derivation of Buchalter et. al. and generalize it to the trispectrum. This applies as 
well to higher-order statistics. The error introduced by the Limber approximation is inferior to 1% for effects with 

i > 200 jHalH. 

Suppose we observe the projection along the line of sight of a three-dimensional statistically homogeneous and 
isotropic random field / 

/>oo 

p{6) = / dmWivfo (Al) 
Jo 

where 9 = (9i, #2,0) and 9 = (9i,9 2 ,l). We propose ourselves to find a relation between the spatial bispectrum 
B p (Ki, K 2 , K3) of p and the spatial bispectrum Bf(K\, K 2 , K 3 ) of /. Following Kaiser method, we consider p to be 
the sum of the contributions from narrow shells with a width Arj much bigger than the relevant wavelength, that is 
9 << At] /t] << 1. This choice allows us to look at fluctuations on scales much less than the characteristic scale over 
which q varies and to assume that contributions from different shells are statistically independent. We then calculate 
the contributions to the bispectrum from each of the shells. At the end we can sum the power for all the shells relying 
on their statistical independence. 

Assuming that q varies very little along the shell and that the section of the shell is plane-symmetric the contribution 
from the shell of width Aij centered at 770 is 

/ .r )0 +A J) /2 

Ap(9)=q(r l0 ) d»//Mi, V0O2, v) (A2) 

Jrio-An/2 

Decomposing the fields in Fourier space, the spectrum of Ap is 

Ap{K) = q(no) / -£4?f(k) / d 2 9e l ^ ok - K ^ e / d V e lk ^ (A3) 



(27T) 



tjo-Atj/2 



where K — (K x , K y , 0) and k = (k x , k y , k z ) = (k, k z ). We can perform the d 2 9 integral which is {2tt) 2 / rj^&^k — K/ijo) 
and the time integral which is Ar]j (k z Ar)/2) to yield 



Ap(K) = I S/(^.^*-)io(W2) (A4) 



A ?7<?(?7o) I dk z ~ ( K X K y 
Vo J (2tt) 770 ' Vo 

Using < /(ki)/(k 2 )/(k 3 ) >= (2?r) 3 B/(A;i, k 2 , &3)<5f,(ki + k 2 + k 3 ) the three-point spectrum is 



< ApiK^ApiK^ApiKs) > = A??3g 6 3fa) / dk lz f dk 2z f dk 3z B f J^ + kl,J^+kl,J^I- + klj 

Vo J J J V Vo V Vo \ Vo 
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'10 




jo(h z Ar)/2)j (k 2z Ari/2)j {{-k u - k 2z )Ar ] /2)5 2 D {K l + K 2 + K 3 ) (A5) 

The important simplification comes from the fact that the major contribution from the first two Bessel functions comes 
from ki z < I/A77 and k 2z < I/A77. But by assumption Ki/r]o » I/A77, K 2 /rj » l/Ar) and K 3 /?] >> I/A77. 
Therefore, k lz « Ki/rj , k 2z « K 2 /rj and k\ z + k 2z « K 3 /rj . So we can neglect all Fourier modes parallel to 



the line of sight. To a very good approximation B f (^ + k 2 z ,^ + k 2 zl ^ + (k lz + k 2z f) ~ %(§-, & 



K 2 , o K 2 , o K 2 

The integration of the Bessel functions gives 

du \ dvjo(u)j (v)jo(u + v) = it 2 (A6) 



We obtain finally 
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< Ap(if 1 )Ap(X2)Ap(K 3 ) > = 4tt 2 ^^ B/(— , — , — + ^2 + #3) (A7) 
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Summing over the shells, and using < p(K 1 )p{K 2 )p{Kz) >= (2ir) 2 B p (K 1 , K 2 , K 3 )5 2 D (Ki + K 2 + K 3 ) 

B P (K U K 2 ,K 3 ) = fd^B^,^) (A8) 
J ?f 7] r) r\ 

The exact same reasoning can be applied to the calculation of the trispectrum. This time we obtain for the same 
projection 

T p ( Kl ,K 2 ,K 3 ,K 4 ) = (A9) 
J rf V V V V 
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